\(~\)

This notebook helps visualize and explore surface and groundwater observations stored in the National Water Information System (NWIS). It was written to perform automated pulls and interactive mapping and plotting of data as it was being collected and analyzed for a specific project, but it can be modified for any region where there is data in NWIS.

\(~\)

1. Map sites

# watershed screening sites 
screensites <- whatNWISsites(bBox = c(-70.448, 41.626, -70.392, 41.677))
nums <- "504|505|506|507"
screensites <- screensites %>%
  filter(grepl(nums, station_nm))

# local demonstration sites
sites <- whatNWISsites(bBox = c(-70.405, 41.669, -70.395, 41.677))

# site info
baseurl1 <- "https://nwis.waterdata.usgs.gov/usa/nwis/qwdata/?site_no="
baseurl2 <-  "&agency_cd=USGS"

siteinfo <- readNWISsite(sites$site_no) %>%
  dplyr::select(site_no, station_nm, site_tp_cd, dec_lat_va, dec_long_va, dec_coord_datum_cd, alt_va, alt_datum_cd, well_depth_va)%>%
  mutate(id = substr(station_nm, 9, 11)) %>%
  mutate(weblink = paste(baseurl1,site_no,baseurl2,sep = " "))

# reclassify site types for mapping
siteinfo$site_tp_cd[grep("M01", siteinfo$station_nm)] <- "multilevel sampler (MLS)"
clust <- "505-0|517-0|524-0|525-0|526-0"
siteinfo$site_tp_cd[grepl(clust, siteinfo$station_nm)] <- "well cluster"
siteinfo$site_tp_cd[which(siteinfo$site_tp_cd=="GW")] <- "water table well"
siteinfo$site_tp_cd[which(siteinfo$site_tp_cd=="LK")] <- "surface water (pond)"

# add a treatment variable
siteinfo$treatment <- NA
controls <- c(505,508,514,515,520,521)
siteinfo$treatment[which(siteinfo$id %in% controls)] <- "control"
impacts <- c(517,522,523,524,525,526)
siteinfo$treatment[which(siteinfo$id %in% impacts)] <- "impact"
siteinfo$treatment[siteinfo$id == 509] <- "other"
# map all sites by type
pal <- colorFactor("Accent", domain = siteinfo$site_tp_cd)

leaflet(siteinfo, options = leafletOptions(zoomControl = FALSE)) %>%
  addProviderTiles('Esri.WorldTopoMap')%>%
  addCircleMarkers(lng = ~dec_long_va, lat = ~dec_lat_va, weight = 1, color = "gray", radius = 6, opacity = 1, 
                 fillOpacity = 1, fillColor = ~pal(site_tp_cd),
                 popup = ~paste(station_nm,"<br>","<a href='", weblink, "' target='_blank'>", "Link to data</a>")) %>%
  addScaleBar("bottomright") %>%
  addLegend(pal = pal, title = "Site type", opacity = 1, values = ~site_tp_cd, position = "topleft")
# map all sites by treatment
pal <- colorFactor("Accent", domain = siteinfo$treatment)

leaflet(siteinfo, options = leafletOptions(zoomControl = FALSE)) %>%
  addProviderTiles('Esri.WorldTopoMap')%>%
  addCircleMarkers(lng = ~dec_long_va, lat = ~dec_lat_va, weight = 1, color = "gray", radius = 6, opacity = 1, 
                 fillOpacity = 1, fillColor = ~pal(treatment),
                 popup = ~paste(station_nm,"<br>","<a href='", weblink, "' target='_blank'>", "Link to data</a>")) %>%
  addScaleBar("bottomright") %>%
  addLegend(pal = pal, title = "Treatment", opacity = 1, values = ~treatment, position = "topleft")

\(~\)

2. Specify water quality parameters for data retrieval

# parameters     (nutrients samples are all filtered)
temp <- "00010"  # temperature                 (degC)
sc <- "00095"    # SpC                         (uS/cm)
ph <- "00400"    # pH
do <- "00300"    # DO                          (mg/L)
nh34 <- "00608"  # ammonia and ammonium, as N                (mg/L)
no2 <- "00613"   # nitrite as N                (mg/L)
no3 <- "00618"   # nitrate as N                (mg/L)
no23 <- "00631"  # nitrate + nitrite as N      (mg/L)
tn <- "62854"    # total nitrogen              (mg/L)
po4 <- "00671"   # orthophosphate as P         (mg/L)
delH <- "82082"  # delta H2/H1                (per mil)
delO <- "82085"  # delta O18/O16              (per mil)

pset <- c(temp, sc, ph, do, nh34, no2, no3, no23, tn, po4, delH, delO)

3. Retrieve and summarize data

# RETRIEVAL

# pull water quality data
## option A (to be retired)
data <- readNWISqw(siteNumber = siteinfo$site_no, parameterCd = pset)
b <- data %>% dplyr::select(c(site_no, sample_dt, sample_tm, parm_cd, result_va))

# ## option B (suggested going forward)
# data <- readWQPqw(siteNumber = paste0("USGS-", siteinfo$site_no), parameterCd = pset)
# # select columns and rename
# b <- data %>% dplyr::select(c(MonitoringLocationIdentifier, ActivityStartDate, USGSPCode, CharacteristicName, ResultMeasure.MeasureUnitCode, ResultMeasureValue))
# names(b) <- c("site_no", "sample_dt", "parm_cd", "charname", "units", "result_va")
# b$site_no <- sub("USGS-", "", b$site_no)

# combine site info and qw data
a <- siteinfo 
a$station_nm[grep("SHUBAEL", a$station_nm)] <- "SHUBAEL POND"
alldata <- right_join(a,b) %>%
  rename(date = sample_dt) %>%
  mutate(param = parm_cd, yearmon = format(date, "%Y-%m")) 

# create param names from codes 
alldata$param <- dplyr::recode(alldata$param, `00010` = "temp", `00095` = "SpC", `00400` = "pH", `00300` = "DO", `00608` = "NH34", `00613` = "NO2", `00618` = "NO3", `00631` = "NO23", `62854` = "TN", `00671` = "PO4", `82082` = "delH", `82085` = "delO")

# pull water levels data
gwl <- readNWISgwl(siteNumbers = siteinfo$site_no) 

# join site info
dtw <- gwl %>%
  filter(!is.na(lev_va))
dtw <- right_join(siteinfo %>% dplyr::select(station_nm, site_no, id, dec_lat_va, dec_long_va, alt_va, well_depth_va, weblink), dtw)%>%
  rename(date = lev_dt)
# SUMMARIES

# number of sample depths at each location
nwells <- alldata %>% 
  group_by(id) %>%
  summarise(depths = length(unique(well_depth_va))) 

# number of depths and samples at each location, by parameter
allsum <- alldata %>% 
  group_by(id, param) %>%
  summarise(count = n()) %>%
  pivot_wider(id_cols = id, names_from = param, values_from = count)

allsum <- left_join(allsum, nwells) %>%
  dplyr::select(id, depths, temp, SpC, pH, DO, PO4, NH34:NO3, TN) %>%
  arrange(depths)

# estimate depth below water table (satdepth) for each sample  
alldata <- left_join(alldata, dtw %>% select(c(date, site_no, lev_va)))        # first, join depth to water
mls <- alldata %>% filter(site_tp_cd == "multilevel sampler (MLS)") %>%  # for mls, assign depth of most shallow sampled port, -0.5ft, as estimate of depth to water
  group_by(id) %>%
  summarise(dtw = min(well_depth_va)-0.5)
for(i in 1:nrow(mls)){alldata$lev_va[alldata$id==mls$id[i]] = mls$dtw[i]}
alldata <- alldata %>%
  mutate(satdepth = well_depth_va - lev_va)

# subset data for a specific sampling event
mapdate <- "2021-11-01"
zt <- alldata %>%
  dplyr::filter(date > mapdate)

# number of depths and samples at each location, by parameter
ztsum <- zt %>% 
  group_by(id, param) %>%
  summarise(count = n()) %>%
  pivot_wider(id_cols = id, names_from = param, values_from = count)

ztsum <- left_join(ztsum, nwells) %>%
  dplyr::select(id, depths, temp, SpC, pH, DO, PO4, NH34:NO3, TN) %>%
  arrange(depths)

# means at each location, single time
zmean <- zt %>%
  filter(!station_nm=="MA-A1W  517-0036") %>%  # remove duplicate shallow well at A517
  #filter(satdepth < 10) %>%                   
  group_by(id, parm_cd) %>%
  summarise(station_nm=first(id), dec_lat_va=median(dec_lat_va), dec_long_va=median(dec_long_va), weblink=first(weblink), well_depth_va=mean(well_depth_va, na.rm=TRUE), result_va=round(mean(result_va, na.rm=TRUE), 1))

# means at each location, time series
tsmean <- alldata %>%
  filter(!station_nm=="MA-A1W  517-0036") %>%  # remove duplicate shallow well at A517
  #filter(satdepth < 10) %>%
  group_by(id, yearmon, param) %>%
  summarise(station_nm=first(id), dec_lat_va=median(dec_lat_va), dec_long_va=median(dec_long_va), weblink=first(weblink), well_depth_va=mean(well_depth_va, na.rm=TRUE), result_va=round(mean(result_va, na.rm=TRUE), 1), treatment=first(treatment))

# max over sample depths at each location
zmax <- zt %>%
  group_by(id, parm_cd) %>%
  summarise(station_nm=first(id), dec_lat_va=median(dec_lat_va), dec_long_va=median(dec_long_va), weblink=first(weblink), well_depth_va=mean(well_depth_va, na.rm=TRUE), result_va=max(result_va, na.rm=TRUE))

# calculate percentages of N in different forms
ftn <- pivot_wider(zt, id_cols = c(id, site_no:dec_long_va, well_depth_va), names_from = param, values_from = result_va)%>%
  filter(!is.na(TN)) %>%
  mutate(percent_organicN = round(100*(1 - ((NH34+NO23)/TN)), 1), percent_NO23 = round(100*NO23/TN, 1))

kable(ftn %>% arrange(id) %>% dplyr::select(station_nm, site_tp_cd, temp, SpC, pH, DO, PO4, NH34, NO3, NO23, TN, percent_organicN, percent_NO23), digits = 2)
station_nm site_tp_cd temp SpC pH DO PO4 NH34 NO3 NO23 TN percent_organicN percent_NO23
MA-A1W 505-0036 well cluster 12.2 143 5.1 9.30 NA 0.02 3.18 3.18 3.23 0.9 98.5
MA-A1W 508-M01-06WT multilevel sampler (MLS) 13.7 105 5.0 3.00 NA 0.02 3.28 3.28 3.38 2.4 97.0
MA-A1W 509-0019 water table well 16.8 105 6.3 0.50 0 0.02 0.27 0.28 0.35 14.0 80.3
MA-A1W 514-0028 water table well 11.7 138 4.7 5.50 0 0.02 6.14 6.14 6.30 2.2 97.5
MA-A1W 515-0036 water table well 11.2 214 4.9 7.80 NA 0.02 11.90 11.90 12.40 3.9 96.0
MA-A1W 517-0037 well cluster 11.3 134 4.8 5.50 NA 0.15 3.65 3.65 3.97 4.3 91.9
MA-A1W 520-0048 water table well 10.7 134 5.3 7.80 NA 0.02 6.56 6.56 6.86 4.1 95.6
MA-A1W 521-0041 water table well 10.7 133 5.2 9.60 NA 0.02 1.58 1.58 1.68 4.8 94.0
MA-A1W 522-M01-06WT multilevel sampler (MLS) 16.4 65 5.5 0.23 0 0.02 0.18 0.18 0.23 12.6 78.7
MA-A1W 523-M01-06WT multilevel sampler (MLS) 15.4 260 4.5 0.09 0 1.73 14.10 14.10 16.70 5.2 84.4
MA-A1W 524-0039 well cluster 10.5 151 5.0 3.10 NA 0.02 3.05 3.05 3.24 5.2 94.1
MA-A1W 525-0042 well cluster 11.7 183 4.8 9.10 NA 0.02 8.50 8.50 8.66 1.6 98.2
MA-A1W 526-0034 well cluster 11.9 133 5.8 0.67 NA 0.02 0.42 0.43 0.53 15.8 80.4

\(~\)

\(~\)

4. Map water quality parameters and depth to water

\(~\)

# function to make the maps
mapparam <- function(data, palette, title){
  leaflet(data, options = leafletOptions(zoomControl = FALSE)) %>%
  addProviderTiles('Esri.WorldTopoMap')%>%
  addCircleMarkers(lng = ~dec_long_va, lat = ~dec_lat_va, weight = 1, color = "gray", radius = 6, opacity = 1, 
                   fillOpacity = 1, fillColor = ~pal(result_va),
                   popup = ~paste(station_nm,"<br>","Value",result_va,"<br>","<a href='"
                                  , weblink
                                  , "' target='_blank'>"
                                  , "Link to data</a>")) %>%
  addScaleBar("bottomright") %>%
  addLegend(pal = palette, title = title, opacity = 1, values = ~result_va, position = "topleft")
}

4a. Mean parameter values

z <- zmean

ttl <- "Temp (°C)"
mapdata <- z %>% filter(parm_cd == temp)
pal <- colorNumeric("viridis", domain = mapdata$result_va)
m1 <- mapparam(mapdata, pal, ttl)
#m1

ttl <- "SpC (uS/cm)"
mapdata <- z %>% filter(parm_cd == sc)
pal <- colorBin("viridis", domain = mapdata$result_va)
m2 <- mapparam(mapdata, pal, ttl)
#m2

ttl = "pH"
mapdata <- z %>% filter(parm_cd == ph)
pal <- colorBin("viridis", domain = mapdata$result_va)
m3 <- mapparam(mapdata, pal, ttl)
#m3

ttl = "DO (mg/L)"
mapdata <- z %>% filter(parm_cd == do)
pal <- colorBin("viridis", domain = mapdata$result_va)
m4 <- mapparam(mapdata, pal, ttl)
#m4

ttl = "Nitrate, as N <br> (mg/L)"
mapdata <- z %>% filter(parm_cd == no3)
pal <- colorBin("viridis", domain = mapdata$result_va)     
m5 <- mapparam(mapdata, pal, ttl)
#m5

ttl = "NH34, as N <br> (mg/L)"
mapdata <- z %>% filter(parm_cd == nh34)
pal <- colorBin("viridis", domain = mapdata$result_va)
m6 <- mapparam(mapdata, pal, ttl)
#m6

ttl = "Total N (mg/L)"
mapdata <- z %>% filter(parm_cd == tn)
pal <- colorBin("viridis",  domain = mapdata$result_va)   
m7 <- mapparam(mapdata, pal, ttl)
#m7

ttl = "Phosphate, as P <br> (mg/L)"
mapdata <- z %>% filter(parm_cd == po4)
pal <- colorBin("viridis", domain = mapdata$result_va)    
m8 <- mapparam(mapdata, pal, ttl)
#m8

# map depth to water
pal <- colorNumeric("viridis", domain = dtw$lev_va)
m9 <- leaflet(dtw %>% filter(date > mapdate), options = leafletOptions(zoomControl = FALSE)) %>%
  addProviderTiles('Esri.WorldTopoMap')%>%
  addCircleMarkers(lng = ~dec_long_va, lat = ~dec_lat_va, weight = 1, color = "gray", radius = 6, opacity = 1, 
                 fillOpacity = 1, fillColor = ~pal(lev_va),
                 popup = ~paste(station_nm,"<br>","Value",lev_va,"<br>","<a href='", weblink, "' target='_blank'>", "Link to data</a>")) %>%
  addScaleBar("bottomright") %>%
  addLegend(pal = pal, title = "Depth to water", opacity = 1, values = ~lev_va, position = "topleft")

sync(m1, m2, m3, m4, m5, m6, m7, m8, m9, ncol = 3, sync.cursor = FALSE)

\(~\)


\(~\)

4b. Maximum parameter values

\(~\)

z <- zmax

ttl <- "Temp (°C)"
mapdata <- z %>% filter(parm_cd == temp)
pal <- colorNumeric("viridis", domain = mapdata$result_va)
m1 <- mapparam(mapdata, pal, ttl)
#m1

ttl <- "SpC (uS/cm)"
mapdata <- z %>% filter(parm_cd == sc)
pal <- colorBin("viridis", domain = mapdata$result_va)
m2 <- mapparam(mapdata, pal, ttl)
#m2

ttl = "pH"
mapdata <- z %>% filter(parm_cd == ph)
pal <- colorBin("viridis", domain = mapdata$result_va)
m3 <- mapparam(mapdata, pal, ttl)
#m3

ttl = "DO (mg/L)"
mapdata <- z %>% filter(parm_cd == do)
pal <- colorBin("viridis", domain = mapdata$result_va)
m4 <- mapparam(mapdata, pal, ttl)
#m4

ttl = "Nitrate, as N <br> (mg/L)"
mapdata <- z %>% filter(parm_cd == no3)
pal <- colorBin("viridis", domain = mapdata$result_va)     
m5 <- mapparam(mapdata, pal, ttl)
#m5

ttl = "NH34, as N <br> (mg/L)"
mapdata <- z %>% filter(parm_cd == nh34)
pal <- colorBin("viridis", domain = mapdata$result_va)
m6 <- mapparam(mapdata, pal, ttl)
#m6

ttl = "Total N (mg/L)"
mapdata <- z %>% filter(parm_cd == tn)
pal <- colorBin("viridis",  domain = mapdata$result_va)   
m7 <- mapparam(mapdata, pal, ttl)
#m7

ttl = "Phosphate, as P <br> (mg/L)"
mapdata <- z %>% filter(parm_cd == po4)
pal <- colorBin("viridis", domain = mapdata$result_va)    
m8 <- mapparam(mapdata, pal, ttl)
#m8


sync(m1, m2, m3, m4, m5, m6, m7, m8, m9, ncol = 3, sync.cursor = FALSE)

\(~\)


\(~\)

5. Plot time series

\(~\)

ttl = "Nitrate, as N, (mg/L) in groundwater wells and pond"
p <- ggplot(alldata %>% filter(site_tp_cd %in% c("well cluster", "water table well", "surface water (pond)")) %>% filter(parm_cd==no3), aes(x = date, y = result_va, group = station_nm)) +
  geom_point(aes(color=id)) +
  geom_line(aes(color=id)) +
  theme_bw()+ 
  ggtitle(ttl) 

ggplotly(p)
ttl = "Nitrate, as N, (mg/L) in multilevel groundwater samplers (MLS)"
p <- ggplot(alldata %>% filter(site_tp_cd %in% c("multilevel sampler (MLS)")) %>% filter(parm_cd==no3), aes(x = date, y = result_va, group = station_nm)) +
  geom_point(aes(color=id)) +
  geom_line(aes(color=id)) +
  theme_bw()+ 
  ggtitle(ttl) 

ggplotly(p)
ttl = "Nitrate, as N, (mg/L) in wells and MLS ports"
p <- ggplot(alldata %>% filter(param=="NO3" & !is.na(treatment)), aes(x = date, y = result_va, group = station_nm)) +
  geom_point(aes(color = treatment)) +
  geom_line(aes(color = treatment)) +
  theme_bw()+ 
  ggtitle(ttl) 

ggplotly(p)
ttl = "Mean nitrate, as N, (mg/L) at site"
p <- ggplot(tsmean %>% filter(param=="NO3" & !is.na(treatment)), aes(x = yearmon, y = result_va, group = station_nm)) +
  geom_point(aes(color = treatment)) +
  geom_line(aes(color = treatment)) +
  theme_bw()+ 
  ggtitle(ttl) 

ggplotly(p)

\(~\)

6. Plot vertical profiles

\(~\)

profdate <- "2021-06-01"

parameters <- c( "NO3")  #"DO", "pH",  "delO", "temp",  "SpC"

profiledata <- alldata %>% 
  filter(id %in% c(505, 508, 517, 522, 523, 524, 525, 526)) %>%
  filter(param %in% parameters) %>%
  filter(date > profdate) %>%
  filter(!station_nm=="MA-A1W  517-0036") %>%  # remove duplicate shallow well at A517
  #mutate(date = as.character(date)) %>%
  #mutate(result_va = ifelse(param == "SpC", result_va/10, result_va)) %>%
  #mutate(param = replace(param, param == "SpC", "SpC/10")) %>%
  arrange(id)
           
xr <- c(0, 50)
yr <- c(80, 0)

# function to make vertical profiles
zprofile <- function(data, title, xrange, yrange){
  p <- ggplot(data %>% group_by(yearmon), aes(x = result_va, y = well_depth_va)) +
    geom_path(aes(color = yearmon, linetype = yearmon))+
    scale_x_continuous()+ #limits = xrange
    theme_bw()+ 
    scale_y_reverse()+ #limits = yrange
    ggtitle(title) 
  ggplotly(p) 
}

# iterated over sites
p <- htmltools::tagList()

for (i in unique(profiledata$id)){
  
  a <- filter(profiledata, id == i) 
  if(i==508){a <- filter(a, site_tp_cd == "multilevel sampler (MLS)")}
  ttl <- paste0(a$site_tp_cd[1], " at ", i)
  p[[i]] <- zprofile(a, ttl, xr, yr)
  
}
p